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ABSTRACT 

We describe a new algorithm for the "perfect" extraction of one-dimensional spectra from two- 
dimensional (2D) digital images of optical fiber spectrographs, based on accurate 2D forward modeling 
of the raw pixel data. The algorithm is correct for arbitrarily complicated 2D point-spread functions 
(PSFs), as compared to the traditional optimal extraction algorithm, which is only correct for a limited 
class of separable PSFs. The algorithm results in statistically independent extracted samples in the 
ID spectrum, and preserves the full native resolution of the 2D spectrograph without degradation. 
Both the statistical errors and the ID resolution of the extracted spectrum are accurately determined, 
allowing a correct \ 2 comparison of any model spectrum with the data. Using a model PSF similar 
to that found in the red channel of the Sloan Digital Sky Survey spectrograph, we compare the 
performance of our algorithm to that of cross-section based optimal extraction, and also demonstrate 
that our method allows coaddition and foreground estimation to be carried out as an integral part 
of the extraction step. This work demonstrates the feasibility of current- and next-generation multi- 
fiber spectrographs for faint galaxy surveys even in the presence of strong night-sky foregrounds. 
We describe the handling of subtleties arising from fiber-to-fiber crosstalk, discuss some of the likely 
challenges in deploying our method to the analysis of a full-scale survey, and note that our algorithm 
could be generalized into an optimal method for the rectification and combination of astronomical 
imaging data. 

Subject headings: methods: data analysis — techniques: spectroscopic — surveys 



1. INTRODUCTION 

Optical fibers offer a compelling advantage for astro- 
nomical survey spectroscopy. By allowing the light from 
targets over a wide field of view on the sky to be rear- 
ranged into a compact format and fed to any number of 
spectrographs in parallel, they can provide a vast mul- 
tiplex advantage over imaging spectrographs. For this 
reason, fiber technology has been adopted for a num- 
ber of major survey programs such as the Las Cam - 



panas Redshift Survey (LCRS: iShectman et al.l 
the Two-degree Field Survey f2dF: IColles~ ct al 



the Sloan Digital Sky Survey (SPSS: lYork et al 



1996), 



2001), 



2000), 



and the recently init iated Baryon Oscilla tion Spectro- 
scopic Survey rBOSS: ISchlegel et al.ll2(ffl9ah of the SDSS3 
project. The use of fiber spectrographs for faint-object 
spectroscopy has however been restrained by concerns 
over throughput and systematic limitations on the qual- 
ity of subtraction of night-sky emission foregrounds. The 
spectroscopic extractions of the SDSS multi-fiber in- 
strument have established a high standard, but signif- 
icant systematic shortcomings remain. Methods to char- 
acterize and partially remove sky-subtraction residuals 
in SDSS fiber spectra can mitigate the problem s ome- 
what fe.g. lBolton et al.ll2004t IWild k. Hewettj|2005ft . but 
do not substitute for a formally correct sky-subtraction 
model. The faintest spectroscopic galaxy surveys have 
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tende d to make use of multi-slit imaging spectrographs 
(e.g [Cowie et al.|[l99l iSteidel et all 120031: iDavis et all 
120031: 1 Abraham et al.ll2004h iLe Fevre et al.ll2005f ) along 
with sky-subtraction te c hniques such as nod- and- shuffle 
dCuillandre et al.l Il994h iGlazebrook & Bland-Hawthornl 
|2001|) or B-spline-based modeli ng of the two -dimensional 
sky spectra in the slits (e.g., iKelsonl [2003D . Nod-and- 
shufflc techniques in particular are ill-suited to most fiber 
spectrographs, require at least double the detector area, 
and furthermore reduce the background-limited signal- 
to-noise by a factor of 1/V2 due to their subtraction of 
data from data. 

This paper outlines an algorithmic framework for the 
modeling and extraction of optical and near-infrared as- 
tronomical fiber spectroscopy to the limit of photon noise 
and native instrumental resolution. By comparing our 
method with the standard techniques of optimal extrac- 
tion currently in wide use, we identify and resolve key 
systematic barriers to "perfect" extraction. As a re- 
sult, multifiber spectroscopy emerges as a clear and com- 
pelling technique for current and future generations of 
faint-galaxy spectroscopic surveys, even well below the 
brightness of the night sky at all wavelengths. This al- 
gorithm will deliver significant benefits to the reanaly- 
sis of the original SDSS (hereafter SDSS1) archive and 
to the ongoing analysis of the BOSS survey, both for 
core redshift-survey goals and for projects that aim to 
select rare emission-line objects from within the regions 
of the spectrum d ominated by OH line emission from 
the night sky (e.g [Bolton et al.ll20Q4t iWillis et alJl2005b 
iBolton et al1 [2006: Will is et al.ll2006l : lBolton et all 2008). 

In considering this subject, we will make a clear dis- 
tinction between the problems of "calibration" and "ex- 
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traction" : calibration is the description of the way in 
which any set of astronomical and environmental stim- 
uli translate into the responses of the digital detec- 
tors (which we assume here to be pixellated charge- 
coupled devices, or CCDs); extraction is the reconstruc- 
tion of particular stimuli from particular responses. More 
strictly speaking, we view an "extracted spectrum" not 
as a model for the flux of the observed source itself, but 
rather as a properly calibrated one-dimensional compres- 
sion of the instrumental response to an observation of 
that source. When executed and reported correctly, an 
extracted spectrum permits a statistically valid x 2 test 
of an input model spectrum against the extracted pixels. 
This paper specifically illustrates a method for carrying 
out this sort of perfect extraction assuming that perfect 
calibration is available. We do not mean to trivialize the 
problem of calibration, and in our concluding remarks 
we will discuss the relationship of our extraction method 
to current and future spectroscopic calibration regimes. 

The paper is organized as follows. Section [2] frames 
the problem of extraction in terms of image modeling, 
lays out the first part of our algorithm, and compares its 
performance with that of standard extraction techniques 
in terms of the quality of their respective models to the 
raw 2D data. Section [3] addresses the issue of resolution 
and covariance in the extracted spectra, and presents 
the second part of our method, which establishes opti- 
mal properties in both these regards. Section |4] presents 
a more realistic demonstration of our algorithm on sim- 
ulated data, illustrating several additional strengths and 
subtleties of the method. Finally, Section \5\ provides a 
discussion and conclusions. 

We will observe the following conventions in this paper. 
Without loss of generality, we will assume that spectro- 
scopic traces run roughly parallel to CCD columns (i.e., 
"vertically"), with wavelength increasing with row num- 
ber and cross-sectional position increasing with column 
number. We will denote vectors in lowercase bold-face 
type (f) and matrices in uppercase bold-face type (A). 
We assume all errors have a Gaussian distribution, and 
we assume no formal priors on fitted model parameters. 

2. EXTRACTION AS IMAGE MODELING 

The current standard of quality for the extrac- 
tion of optical fiber spectroscopic data from digital 
images to ID spectra is the optimal te chniq ue de- 
scribed by j Hewett et all l| 19851 ). iHornd (|1986[ ). and 
IRobertsonI ( 19861 ). a nd subsequently implement ed 



countless forms (e.g.. iMarsh 



Hall et al.|[l99l iBacon et al 



1989: 



2001 



Kinnev et all fl 991; 
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20021: iCushing et al l "l200l iBershadv et alj 120051; 
Zanichelli et al.|l2005l:lDixon et al.ll2007HBolton fc Buries! 
2007tlCui et al.ll2008f ). This algorithm models the two- 
dimensional spectrum image row by row with a 
minimum-x 2 scaling of the cross-sectional profile. The 
fitted amplitude is taken as an optimal estimator of the 
ID spectrum at the wavelength corresponding to that 
row. Optimal extraction as currently understood and 
practiced has significant advantages, most notably in 
the increased signal-to-noise ratio of the extracted spec- 
trum in comparison with the boxcar-aperture summing 
technique. It also allows for model-based flagging of 
pixels afflicted by cosmic-ray hits. Optimal extraction 
furthermore leads to statistically uncorrelated extracted 



spectrum values, although this property is not preserved 
under the rebinning that is required for the combination 
of multiple exposures that are not precisely aligned at 
the raw-pixel level. 

The most significant (and under-appreciated) short- 
coming of the optimal extraction method is that it is only 
correct in the case where the two-dimensional image I of 
a monochromatically illuminated fiber (the point-spread 
function, or PSF) is a separable function of column x 
and row y offset from the PSF centroid: 



I(x,y) = I x (x)I y (y) 



(1) 



If this condition does not hold — and it does not hold for 
nearly all PSF models beyond a single Gaussian ellipsoid 
whose principal axes are aligned with the CCD rows and 
columns — then the cross-sectional profile of the 2D spec- 
trum is no longer fixed, but rather depends upon the 
details of the input spectrum itself. This is of limited 
concern if the spectrum is smoothly varying with wave- 
length. However if the spectrum involves many sharp fea- 
tures such as the bright OH rotational emission lines in 
the near infrared spectrum of the night sky, then the out- 
put of the row-by-row optimal extraction algorithm be- 
comes ill-defined. For example, the cross section through 
the core of an emission line may have a different shape 
than the cross section through the wing of the line, yet 
the extraction will attempt to model them both with a 
single average cross-sectional shape. This will degrade 
the resolution and signal-to-noise ratio of the extracted 
spectrum, and may introduce a bias into the spectrum 
estimation. 

The shortcomings of the row-wise optimal extraction 
method can be overcome by modeling the 2D spectro- 
scopic data image in a manner that more accurately re- 
flects the way in which an input spectrum translates into 
photon counts in the CCD image. For the case of a linear 
CCD detector (which we will assume), the system cali- 
bration can be described as a matrix A, whose elements 
An describe the (noise-free) counts in CCD pixel i for 
a unit of monochromatic input at wavelength A^. Note 
that we are suppressing the natural two-dimensionality 
of the CCD by allowing i to index all column/row combi- 
nations. We are also suppressing any possible multi-fiber 
dependence; it is straightforward to incorporate this sort 
of dependence into the ^-indexed dimension, and we will 
in fact do so in £|4] below. The calibration matrix A in- 
corporates all the effects of wavelength calibration, spec- 
trum trace position, PSF shape and its dependence upon 
position, flat-fielding effects, spectrophotometry, and ex- 
posure time. Note that A will generally be a sparse ma- 
trix, since an input at a given wavelength in a given fiber 
will only influence a relatively small number of the full 
CCD pixels. For the case of an input spectrum vector f , 
the observed CCD pixel count vector p will be given by 



p = Af 

or in indexed notation 



n 



Pi 



(2) 



(3) 



where n is a pixel noise vector. Extraction is then the 
reconstruction of f from p given a knowledge of A. Al- 
though extraction is non-trivial, it is a linear problem. 



Spectro-Perfectionism 



3 



Since the wavelength coordinate is continuous, Equa- 
tion [2] should most properly be written in integral rather 
than matrix form. However our interest is in solving 
for f: inversion for an infinite number of differentially 
spaced wavelength amplitudes from a finite number of 
pixels and resolution elements is neither possible nor de- 
sirable. We therefore restrict our attention to a discrete 
set of wavelength sampling positions {A^}. The means of 
determining the most appropriate sampling density for 
these positions will be addressed in more detail further 
below; in general they may be either more closely or more 
widely spaced than the pixel rows in the raw data. How- 
ever, if they are too closely spaced, then the solution for 
f becomes degenerate. 

Assuming Gaussian noise, the minimum-^ 2 solution for 
the input spectrum vector f from the data vector p is 



f = 



(A T N _1 A) 



A T N _1 p . 



(4) 



Here, N is a pixel noise matrix, with iVy = (niiij). We 
treat raw pixel errors as statistically independent, and 
thus the noise matrix is diagonal. Our model for the 
2D raw pixel data is simply given by m = Af . Again, 
this requires that the sampling points not be too closely 
spaced in wavelength, so that A T N -1 A does not become 
non-invertible. 

To demonstrate how this extraction method improves 
upon row-by-row optimal extraction based on a fixed 
cross-sectional profile, we will adopt the following illus- 
trative fiber spectrograph PSF: 



I{x,y) = 



(1-6) 



exp 



(2<r») 



bexp(-r/r ) 
27rror 



(5) 



Here, a (in pixels) controls the size of a Gaussian core 
to the profile (first term), and tq sets the characteristic 
size of the profile wings (second term). The parame- 
ter b controls the fraction of the total flux in the wing 
component. The form of the wing component is taken 
from a parameterization of the near-IR scattering seen 
within the SDSS1 spectroscopic CCDs (J.E. Gunn, un- 
published). The coordinate r — \Jx 2 +y 2 is the radial 
offset in pixels on the CCD from the center of the PSF 
spot. The form r c ii in the argument of the Gaussian term 
indicates that the core profile can have an elliptical shape 
as determined by 



fell 



qx 



y' 2 /q , 



(6) 



where q is a minor-to-major axis ratio, and (a/, y') are 
related to the CCD column/row coordinates (x,y) by 
a possible rotation and translation. For the wings, we 
adopt the values b — 0.1 and r$ = 36 pixels, roughly ap- 
propriate to the SDSS1 CCDs at a wavelength of 8500 A. 
For the core, we take a = 1.1 pixels, roughly the median 
characteristic value for the SDSS1, and an ellipticity of 
q = 0.75 with a major axis position angle inclined at 
45° relative to the column/row axes. This core ellip- 
ticity is somewhat more exaggerated that what is seen 
in well-focused SDSS1 spectrograph exposures, but will 
provide a good test of our method in the presence of the 
astigmatic aberrations that can arise in wide-field camera 
systems. 

Generally speaking, the calibration matrix A is de- 
termined by the combination of the fiber PSF func- 



tional form, the dependence of its parameters upon wave- 
length and position (which can be derived from calibra- 
tion frames illuminated by gas discharge lamps that have 
multiple discrete emission lines at known wavelengths), 
and the wavelength- and pixel-dependent throughput of 
the system (which can be derived from calibration frames 
uniformly illuminated by incandescent lamps with ap- 
proximately flat continuum spectra). For simplicity, we 
will use CCD pixel rows as a proxy for wavelength, ne- 
glect any variation of the PSF shape parameters with 
position, and neglect any pixel-to-pixel sensitivity varia- 
tions. Assuming that all these effects can be calibrated 
successfully, they can be incorporated into the approach 
that we describe. We also neglect the problem of cos- 
mic rays here, although a model-based extraction such as 
ours is ideally suited to the iterative masking of cosmic- 
ray pixels in the raw data based upon their statistical 
disagreement with the best extraction model. Finally, 
we ignore the possibility of a large-scale scattered-light 
component in the image, but the incorporation of this 
phenomenon (as well as those of imperfectly subtracted 
instrumental bias and dark-current levels) into our 2D 
modeling would be straightforward. 

The most stringent test of any extraction is provided 
by unresolved emission lines, which show the sharpest 
real features possible in the data. To test our method 
alongside the standard row-wise extraction, we therefore 
create a noise-free spectrograph-like image consisting of 
11 emission lines of equal flux spaced vertically over 80 
pixels within a 31x101 pixel image. We include a slight 
tilt of the spectral trace relative to the vertical, and in- 
corporate a mild rate of change in the spacing of the lines 
so as to ensure a range of sub-pixel positions for the line 
centers. Each pixel value is computed by integrating the 
profile with 4x4 sub-sampling. We incorporate a small 
core radius in the denominator of the scattering wing 
term to avoid numerical difficulties with the (integrable) 
singularity. This simulated image is shown in the left- 
most panel of Figure [1] 

Next we perform an extraction of this model image us- 
ing both the standard row-wise, cross-section based opti- 
mal extraction and our 2D modeling method. For the 2D 
extraction, we compute a basis set of 101 PSFs centered 
on the trace position in each of the 101 image rows. Note 
that the "emission lines" of the simulated image will in 
general not be coincident with these basis functions, but 
will instead be offset from the nearest basis PSF by some 
sub-pixel amount in y. For the row-wise extraction, we 
sum all 101 basis functions with equal weight to gener- 
ate the image of a flat continuum spectrum, and take 
the normalized cross-sectional shape of this spectrum in 
each row as the basis profile for the extraction of that 
row. Thus both the 2D and row-wise extractions will 
have an equal number of free extraction parameters, with 
the same sampling interval between them. Giving equal 
weight to all pixels, we compute the linear least-squares 
set of amplitudes that multiply the basis functions of 
each extraction method to best reproduce the simulated 
image. The results are shown in Figure [1] The differ- 
ence between the quality of the two models can be seen 
immediately, with the row-wise extraction model unable 
to capture the variation of the cross sectional shape with 
spectrum that comes from the ellipticity of the PSF core 
and the presence of the scattering wings acting in con- 
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Fig. 1. — Left: from left to right: simulated noise-free emission-line image, 2D extraction model of simulated image, and row-wise 
extraction model of simulated image. Color scaling is in units of the base-10 logarithm of the pixel value, with the pixel values themselves 
scaled to have an average value of unity across the entire image. Right: residuals resulting from the subtraction of the 2d extraction model 
(left) and row-wise extraction model (right) from the simulated emission-line image. Color scale is in (non-logarithmic) units of residual 
counts. Image regions are 31x101 pixels in size in all cases. 



cert with the sharp spectral features. Figure[T]also shows 
the residuals of subtraction of the two extraction mod- 
els from the simulated image. The residuals of the 2D 
method have a scale roughly 200 x smaller than the resid- 
uals of the row-wise method. 

We can make a quantitative comparison of the relative 
quality of the two extraction models to the simulated 
image in terms of the following two figures of merit: 



D 




(7) 



(8) 



find D sq = 0.0006 and D ab 
and -D sq = 0.12 and D a b = 



Here, the pi are the raw pixel data values, and the rrij are 
the 2D extraction model values in those pixels. Both of 
these figures will be equal to zero in the case of a perfect 
model, and equal to 1 for a model that is identically zero 
(i.e., no model at all). For the case described above, we 

= 0.0008 for the 2D model, 
0.20 for the row-wise model, 
again confirming a factor of w200 improvement with the 
2D model. 

Before proceeding further, we address the question 
of appropriate sampling. While the sampling of the 
row-wise extraction is naturally limited to one sampling 
point per row, the sampling of the 2D extraction method 
can easily be adjusted to be either finer or more coarse 
than this. The residuals seen above after subtracting 
the 2D model to the simulated emission-line image, al- 
though very small compared to those of the row-wise 
model, nevertheless exceed the numerical precision of the 
calculation, and are due to the misalignment between 
the extracting basis and the input emission lines in the 
"wavelength" direction. We might reasonably improve 
upon the extraction model by adopting a more finely- 
spaced basis set. If the typical fiber spectrograph PSF 
were purely band-limited, then the appropriate sampling 
could be chosen via the Nyquist criterion. However, fiber 



spectrograph PSFs are primarily determined not by a 
diffraction limit, but rather by the convolution of a sharp 
"top-hat" fiber image with a set of optical aberrations, 
and thus the determination of an appropriately critical 
extraction sampling must be made according to more 
empirical factors. For the particular case under consid- 
eration here, we have tested the quality of the 2D model 
to the simulated image for sampling densities ranging 
from 0.5 to 2 basis functions per row in the 2D image. 
Both D sq and D a b improve with increasingly fine sam- 
pling until about 1.3 basis functions per row, where they 
attain a value of a few xlO -5 . With finer sampling be- 
yond this point, the model becomes worse as the ma- 
trix to be inverted approaches singularity. In practice, 
the appropriate sampling density will be determined by 
the details of the spectrograph PSF and dispersion, and 
consequently by the spectral resolution as defined in the 
following section. 

3. RESOLVING THE RESOLUTION 

We have demonstrated the fidelity of our 2D extrac- 
tion model to the two-dimensional pixel data, but what 
of our extracted spectrum? Naively, f from Equa- 
tion [?] is our ID spectrum estimator, and the matrix 
C = (A T N _1 A) _1 is the covariance matrix of its in- 
dividual values. But all is not well: Equation 0] not only 
extracts the spectrum, it deconvolves the spectral reso- 
lution. The instability of the deconvolution manifests 
itself in a significant "ringing" in the extracted spec- 
trum, along with large correlations and anti-correlations 
between extracted pixels. Both of these are decidedly 
undesirable features in a spectrum. The way to remedy 
this situation is the second key element of our extraction 
method: to re-convolve our de-convolved spectrum back 
to the same resolution as the raw data. 

Consider the inverse covariance matrix C _1 of the de- 
convolved spectrum basis: 

Cr 1 =A T N -1 A. (9) 
This matrix is symmetric and band-diagonal, and all of 
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its entries are non-negative (assuming that there is no 
way for any input spectrum to subtract counts from the 
CCD). Now take the unique non- negative matrix square 
root of C _1 to find a symmetric matrix Q such that 

CT^QQ. (10) 

This may be done by determining the eigenbasis of C _1 , 
taking the element-wise square root of the diagonal ma- 
trix of its eigenvalues (which will all be positive since 
C _1 is positive definite), and transforming this new di- 
agonal matrix back using the unitary matrix that relates 
the eigenbasis to the original basis. Equation [TOl would 
still hold if we applied any arbitrary set of sign-flips to 
the square-root-eigenvalue matrix (leading to other ma- 
trix square roots besides the unique non- negative one), 
but these would introduce undesired negativity into Q. 
Next, define a normalization vector s through 

s e = ^2Qie' ( n ) 

i 

a matrix R through 

R- u = s~ x Q- u (no sum), (12) 

and a diagonal matrix C^ 1 with entries given by 

Cjl = sj . (13) 
By construction, we now have 

CT^R^CH-R, (14) 

and consequently 

C = RCR T . (15) 

Note the complete analogy between Equations [5] 
and 1141 The matrix R defines a transformation from 
the deconvolved spectrum to a "reconvolved" spectrum 
for which (1) the extracted ID pixels are statistically in- 
dependent of one another (since the matrices C _1 and 

C are diagonal), and (2) the blurring of input spectra 
in wavelength from the deconvolved basis is statistically 
equivalent to the blurring imposed by the actual 2D spec- 
trograph matrix A on the physical input spectrum. It 
is in this sense that we can re-convolve with the "same" 
resolution as was inherent in the 2D data. Our extracted 
ID spectrum is then 

f = Rf, (16) 

with uncorrelated errors described by the diagonal co- 
variance matrix C and undegraded resolution described 
by the matrix R. We may think of this reconvolved spec- 
trum as a model for what would have been observed by a 
truly one-dimensional spectrograph with the same reso- 
lution as our two-dimensional CCD spectrograph. Since 
the resolution of the extracted spectrum is precisely char- 
acterized by the matrix R, we can convolve any theoret- 
ical model for the input ID spectrum with this matrix 
in order to compare to the data and compute a statis- 
tically correct x 2 value. The sense of the normalization 
imposed by Equation [12] is to describe the elements of 



f as a weighted sum over the elements of the theoreti- 
cal input spectrum, with the sets of weights individually 
summing to 1. 

Figure [2] illustrates the deconvolved and reconvolved 
spectra that result from the extraction of the simulated 
emission-line spectrum of $21 Note that the ringing in the 
deconvolved spectrum is completely absent from the re- 
convolved spectrum. The peaks in the reconvolved spec- 
trum represent the true undegraded resolution of the 2D 
spectrographic data. The values in the deconvolved spec- 
trum are highly correlated with one another, whereas the 
values in the reconvolved spectrum are completely uncor- 
related. Also shown is a representation of the difference 
between the extracted spectra using the 2D and row-wise 
methods. Although the two spectra are sufficiently simi- 
lar that their differences cannot be seen in a direct plot, 
they differ at the few percent level, which is significant 
for application to the extraction of faint galaxy spectra 
in the presence of strong night-sky line emission. Also 
note that the resolution of the 2D-extracted spectrum is 
higher than the resolution of the row-extracted spectrum. 

In defining the resolution with which to reconvolve our 
deconvolved extracted spectrum, we have explicitly cho- 
sen a form that results in an exactly diagonal sample co- 
variance matrix in the reconvolved basis. Additionally, 
as a consequence of the band-diagonal nature of C^ 1 , our 
R matrix is also essentially band-diagonal, correspond- 
ing to a localized "line-spread function" in the extracted 
ID spectrum as seen in Figure [2j It is nevertheless worth 
remembering that other choices for R are available. For 
instance, there may be applications for which it is desir- 
able to have a simpler parameterized form for the res- 
olution than will in general result from the procedure 
outlined above. The parameters of this form could be 
optimized so as to make the off-diagonal correlations in 
the reconvolved basis sufficiently small while still non- 
zero. 

4. A MORE INTERESTING TEST 

Before taking up the longer-term challenge of imple- 
menting our algorithm on actual survey spectroscopy, we 
illustrate its power with a more realistic simulation than 
the simple case above. In particular, we now include the 
following effects, all of which will be found in real survey 
data: 

1. Noise, 

2. A varying 2D fiber PSF shape, 

3. Overlap ("cross-talk") between the 2D images of 
neighboring spectra on the CCD detector of a mul- 
tifiber spectrograph, 

4. Multiple exposures with relative movement ("flex- 
ure") of the fiber images on the CCDs between 
them, and 

5. A night-sky foreground that may vary between ex- 
posures and must be modeled and subtracted to 
reveal a much fainter object spectrum of interest. 

We use nearly the same PSF model as above, but now 
generate four "spectra" on the same image, separated by 
5.7 pixels from one another in the horizontal direction. 
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This separation is somewhat smaller than the «6.2 pixel 
separation between neighboring spectra in the SDSS1, 
and therefore leads to greater cross-talk. The Gaussian 
core of each fiber is given a different ellipticity (ranging 
from 0.9 to 0.6) and position angle (ranging between 15° 
and 60° from the x pixel direction) , to mimic optical vari- 
ations and fiber heterogeneities. The central two spectra 
are taken as "object" spectra, and the outside two are 
"sky" spectra. We also generate three "exposures" of 
these four spectra, shifting the fiber images relative to 
the CCD grid in each case to simulate spectrograph flex- 
ure. We generate a "sky" spectrum of 15 emission lines 
at increasing separation over the range of the spectra. In 
each of the three exposures, the amplitudes of these indi- 
vidual sky lines are randomized, and the sky realization 
for each exposure is projected through all four "fibers" in 
that image. The positions of the sky lines in wavelength 
are fixed across all exposures. Next, two "object" spectra 
are projected through the central two fibers of all three 
exposures and added to the image of the sky spectra. 
The relative fluxes are taken so that, in the fibers con- 
taining both object and sky, the total sky counts exceed 
the total object counts by a factor of 20. We use an ex- 
traction sampling density of 1.2 basis functions per CCD 
row. The noise level is set by assigning 10 4 statistical sky 
counts (i.e., photons) per spectrum per sampling point 
in total across all three exposures. This implies 500 total 
object counts per sampling point, and an approximate 
background-limited signal-to-noise ratio of 5 per extrac- 
tion pixel for the objects. We also include a "read-noise" 
term of 5 counts per pixel RMS. We use the statistical 
noise level to generate a noise-image realization that we 
add to the sky and object data. The resulting three sets 
of four spectra are seen in the leftmost third of Figure S 
The full power of the 2D modeling extraction now be- 
comes apparent. We construct a generalized calibration 
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Fig. 2. — Extracted spectra of simulated noise-free emission-line 
image. Thin black line: deconvolved spectrum from the 2D mod- 
eling extraction method. Note significant ringing. Thick blue 
line (upper of the two thick lines, rendered with steps): decon- 
volved spectrum from 2D modeling extraction reconvolvcd to the 
native spectrograph resolution using the resolution matrix R de- 
fined in the text. Thick green line (lower of the two thick lines, 
rendered with steps): 10 X the difference between the reconvolved 
2D-model extracted spectrum and the row-model extracted spec- 
trum. "Upward-peakiness" at the position of the emission lines in- 
dicates that the 2D-extracted spectrum has higher resolution than 
the row-extracted spectrum. 



matrix, involving input spectra of both the night sky (one 
for each exposure) and the science target objects (one for 
each object). In each exposure, the sky basis projects 
into the images of all fibers — sky and object alike — while 
the object bases project only into the images of the object 
fibers. The multiple exposures of each object spectrum 
are extracted in a single step: there is no need for the 
separate steps of registration, extraction, and coaddition 
of individual frames. Similarly with the sky spectra, we 
extract a single sky spectrum from all fibers of a given ex- 
posure. Within the object fibers, we extract the sky and 
object spectra in sum together. There is no explicit sky- 
subtraction step, but rather a sky-object decomposition 
that is an integral part of the extraction into individual 
ID spectrum components. All at once, we model and 
extract three skies — one for each exposure — and two ob- 
ject spectra. In summary: extraction, coaddition, and 
sky subtraction are all subsumed into a single image- 
modeling operation. The results of this modeling can be 
seen in the second and third parts of Figure [3l illustrat- 
ing the sky-object decomposition and the "photon-noise" 
limited quality of the sky+object model to the three 2D 
exposures. 

We note that the all-in-one sky modeling and extrac- 
tion is necessary in order to model and decompose the 
sky spectrum from the object spectra in the deconvolved 
frame, upstream from the fiber-to-fiber PSF variations 
that would lead to systematic errors in a traditional 
modeling and subtraction of the sky from extracted and 
resolution-convolved spectra. If necessary, an accounting 
for any spatial variations of the sky brightness over the 
telescope focal plane could be built in by way of addi- 
tional linear sky components. The combination of multi- 
ple object spectra at the time of extraction, on the other 
hand, is not strictly necessary. Sky-subtracted object 
spectra from the individual exposures could be combined 
together in a subsequent step, to allow for the non-linear 
determination of spectrophotometric variations between 
the frames. However, once these variations are deter- 
mined, the proper combination of these multiple frames 
would resemble a second extraction, with the concatena- 
tion of their individual R matrices forming the new A 
matrix (in the notation of $2J) . 

The presence of crosstalk between neighboring spec- 
tra, and between the object spectra and the skies, makes 
the determination of the ID resolution matrix somewhat 
more challenging. Consider the full covariance matrix C 
and inverse covariance matrix C _1 in the deconvolved 
extracting basis: they consists of multiple blocks, with 
each block describing the statistical coupling between 
different extracted spectra. The matrix C _1 is band 
diagonal in all blocks, but it has non-zero elements in 
the off-diagonal blocks coupling different extracted spec- 
tra to one another. Assuming the spectral cross-talk on 
the CCD is from spectra that are otherwise unassociated 
with one another, we do not simply want to take the 
matrix square root of C _1 to define our resolution, since 
the resolution defined in this way would mix extracted 
spectra with one another. A possible solution which we 
adopt here is to invert C _1 to obtain C, zero out all 
the entries in the off-diagonal blocks of C, re-invert, and 
define the resolution in terms of the square root of the 
resulting matrix. To make this explicit with block matrix 
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Fig. 3. — Simulated multifiber, multi-exposure spectroscopic data, including noise, flexure, a non-uniform PSF, and "sky". Left: Three 
"exposures" of four fiber spectra, including simulated flexure and sky-spectrum variability. Object spectra are included in central two fibers 
in each set, but are too faint to see directly. Center: same as left, after subtraction of extraction model for the sky component in each 
exposure, with gray-scale stretched by a factor of 40 to reveal the traces of the object spectra. Right: As in center, but after subtraction of 
extraction model for the object spectra as well. When scaled by the pixel errors, these residuals are consistent with a reduced x 2 °f unity. 
Each "exposure" is 51x101 pixels in size. 



notation for a case of two extracted spectra: 

-l 



c=(c- 1 



Cn c 12 

C21 C22 



Cu 

c 22 



v- 1 



[C11] 






^22] 



QQ 



(17) 
(18) 
(19) 



and proceed to define the resolution matrix R in terms 
of this Q as indicated by Equations [TT1 and [T2"l Note that 
[Cu] -1 ^= C^f and [C22] -1 ^ C22 due t° the non-zero 
off-diagonal blocks in the constructed inverse covariance 
matrix C _1 . 

It can be shown that by this definition, the resolution 
matrix will diagonalize C within its diagonal sub-blocks 
so that extracted samples will be statistically indepen- 
dent of one another within each spectrum. There will in 
general be non-zero correlations between the extracted 
samples of different spectra with one another, but this is 
unlikely to be important for spectra that are otherwise 
unrelated. Note that if the problem were rather to ex- 
tract the spectra of an integral-field spectrograph, with 
cross-talk between fibers that were also adjacent to one 
another on the sky, then it would perhaps be desirable in- 
stead to have a resolution matrix that mixes spectra with 
one another but provides full statistical independence of 
samples both within a single spectrum and among neigh- 
boring spectra. 

Computing the resolution of our simulated multi-fiber, 
multi-exposure set in the manner outlined above, and 
re-convolving the extracted object and sky spectra, we 
obtain the results shown in Figure |4] The sky spectra 
are scaled down by a factor of 20 for display purposes. 
When weighted by the error estimates, the extracted 
spectra have a reduced x 2 of approximately unity relative 
to the resolution-convolved input spectra. This demon- 
strates the suitability of our approach for the extraction 
of faint galaxy spectra in the presence of high and sharply 
wavelength-dependent night-sky foregrounds. 

5. DISCUSSION AND CONCLUSIONS 

The extraction algorithm that we have described does 
not come without a price. In the presence of fiber-to-fibcr 
cross-talk, the standard row-wise optimal extraction will 
couple extracted amplitudes between fibers in a single 
row, leading to a banded matrix of dimension equal to 
the number of fibers that must be inverted; this process 



must then be repeated for each row in each exposure. 
Our 2D modeling extraction, on the other hand, couples 
extracted amplitudes between fibers, wavelengths, and 
exposures. Thus the matrix to invert for the solution set 
of spectra has sides of dimension equal to 

^spectra x ^samples per spectrum x ^exposures ■ (20) 

For one SDSS1 pointing, this would correspond to 320 
fibers (one of the two spectrographs), approximately 
4000 sampling points, and three exposures: i.e., a square 
matrix nearly 4 million to a side. With a brute-force 
approach this matrix could never be stored, let alone in- 
verted, with any hardware of the present or foreseeable 
future. The way forward to determining the extracted 
spectra will no doubt lie in exploiting the sparseness of 
the inverse covariance matrix to reduce storage and com- 
putation, and to apply an iterative method such as the 
conjugate gradient to solve for the extracted spectra. To 
determine the resolution with which to reconvolve the ex- 
tracted spectrum, which formally requires the inversion 




40 60 
Pixel row 

Fig. 4. — Extracted spectra of two simulated objects from multi- 
ple exposures as described in ij4]and depicted in Figure [3] together 
with extracted "sky" spectra. Black lines (rendered with steps) 
show the extracted spectra, while blue lines (rendered smoothly 
and tracing the black lines) show the input spectra convolved with 
the ID resolution. Green lines of varying shades (thinner, with 
higher peak values, and identical in both plots) indicate the ex- 
tractions of the three different realizations of the "sky" spectrum 
in the three individual exposures, divided by a factor of 20 for 
display purposes. 
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of the full inverse covariance matrix C _1 , the practical 
solution will be to invert a sufficient subspace of C _1 
surrounding each spectrum (or subsegment of a spec- 
trum) to define an acceptably accurate approximation 
to the desired resolution matrix. The exact requirements 
will depend upon the degree of cross-talk between neigh- 
boring fibers in the spectrograph under consideration. 
Even with these strategies, a usable implementation of 
our algorithm for real multifiber survey data may require 
high-performance parallel computing, depending on the 
computational expense of the necessary matrix-element 
calculations. 

Most immediately, we plan to implement the strategy 
outlined in this paper to the extraction of spectra from 
the BOSS Survey. We also plan to conduct a reanalysis 
of the SDSS1 archive, to provide the definitive version of 
this important spectral database with the best possible 
extracted resolution, signal-to-noise, and foreground sub- 
traction. These techniques also offer promise for the up- 
coming Apache Point Observatory Galaxy Evolut ion Ex- 
periment ( APOGEE: lAllende Prieto et all 12003 ) of the 
SDSS3, a high-resolution, near-infrared multifiber sur- 
vey of red giant stars in our Galaxy that aims to con- 
strain their evolutionary history as traced by multiple 
chemical abundances. This extraction strategy will also 
form a key pa rt of the technical fe asibility of the Big- 
BOSS survey (|Schlegel et alj l2009bh . which proposes to 
use a 5000-fiber spectrograph fed by a 3° field-of-view 
focal plane positioner system on a 4m-class telescope 
to measure the baryon acoustic scale and redshift-space 
distortions over the redshift range 0.2 < z < 3.5. Al- 
though the implementations for these different surveys 
will differ in detail, we believe that the software engine 
for the core extraction computations can be written in a 
general-purpose form. The application of this technique 
to slit spectroscopy may be possible as well, although 
the preservation of object spatial information by the slit 
makes the problem substantially more complicated. 

In all cases, the full power of this extraction can only be 
realized with sufficiently accurate calibration. The cur- 
rent standard calibrations for fiber spectroscopy are ex- 
posures of uniform spatial illumination by flat-spectrum 
incandescent lamps ( "flats" ) and multi-emission line gas- 
discharge lamps ("arcs"). Within the framework of our 
extraction algorithm, the former will determine the rel- 
ative sensitivity of the individual fibers and pixels, and 
the latter will determine the spectrograph PSF shape 
and position as a function of illuminating wavelength in 
each fiber. Assuming there are no systematic offsets be- 
tween the calibrations and the science exposures, and 
assuming that the variation of the spectrograph PSF is 
smooth enough with wavelength to be well-sampled by 
the arc frames (which are sparsely distributed in wave- 
length), these calibrations should contain sufficient in- 
formation to determine the calibration matrix A. We 
may proceed by extracting the arcs and flats together, 
with each one described by a single spectrum projected 
through an initial guess for A, and then optimize the 
parameters of A by non-linear iteration so as to improve 
the quality of the 2D extraction models to convergence. 
A more direct and ambitious approach to the determi- 
nation of the calibration matrix would be to illuminate 
the facility calibration screen with eith er a high-wattage 
monochrometer or a tunable laser (c.f. iStubbs fc Tonrvl 



Schlegel 

120061 : iCramer et al.l 12009ft , and to step the illumination 
source through wavelength in subsequent exposures so as 
to map out A explicitly. In practice, with the exception 
of spectrograph systems that are very stable thermally 
and mechanically, there is likely to be some shifting of 
the fiber positions and focus on the CCDs between the 
calibration and science frames. In this case, "tweaks" 
to the parameters of A will be derived from the shapes 
and positions of the fiber traces and night-sky emission 
lines in the science frames. Finally, we note that the cal- 
ibration may be significantly improved by incorporating 
all available knowledge about the optical design of the 
telescope and instrument, rather than treating the sys- 
tem as a black box to be specified entirely by empirical 
calibration data. 

Putting aside the computational challenges that arise 
from the consideration of continuously two-dimensional 
input data, the method we have described can also be 
generalized into an optimal recipe for the rectification 
and combination of multiple CCD imaging exposures: 
i.e., taking the f of £|2] to be a two-dimensional im- 
age model to be extracted from the data. As with 
the spectroscopic application, the resulting extractions 
would have optimal resolution, statistically independent 
extracted image samples, and a definition in terms of the 
optimization of a well-motivated scalar objective func- 
tion describing the quality of a model fit to all of the in- 
dividual exposures. The implementation would be non- 
trivial, but the benefits could be great. A significant 
challenge on the calibration side is that the details of the 
imaging PSF, which must be known, will depend upon 
the spectral energy distributions of the imaged objects, 
which will in general be unknown. 

In summary, we have demonstrated a new spectrum 
extraction algorithm for optical and near-infrared astro- 
nomical fiber spectroscopy. Given sufficiently accurate 
calibration, this method can extract spectra to the sta- 
tistical noise limit in the presence of arbitrarily compli- 
cated point-spread functions and arbitrarily high and 
wavelength- varying foregrounds. The extracted spec- 
tra have optimal and precisely quantified resolution and 
signal-to-noise, along with statistically uncorrelated pix- 
els, for any number of sub-exposures in combination to- 
gether. As such, statistically accurate x 2 tests may be 
made between the extracted data and theoretical mod- 
els of the input object spectrum. This algorithm rep- 
resents a fundamental improvement upon the current 
state-of-the-art methods in use for the extraction of fiber 
spectroscopy, and thus motivates a serious and positive 
reevaluation of the promise of fiber-fed spectrographs for 
next-generation ground-based faint-object surveys. 

The authors wish to thank Scott Buries, Julian Borrill, 
Robert Lupton, David Hogg, Sam Roweis, and Michael 
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subject. DJS acknowledges the support of the Director, 
Office of Science, of the U.S. Department of Energy under 
Contract No. DE-AC02-05CH11231. 
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